arXiv:math/9808109vl [math.NA] 26 Aug 1998 


SPATIAL DISCRETIZATION OF PARTIAL 
DIFFERENTIAL EQUATIONS WITH INTEGRALS 

ROBERT I. MCLACHLAN 


Abstract. We consider the problem of constructing spatial finite 
difference approximations on a fixed, arbitrary grid, which have ana¬ 
logues of any number of integrals of the partial differential equation 
and of some of its symmetries. A basis for the space of of such differ¬ 
ence operators is constructed; most cases of interest involve a single 
such basis element. (The “Arakawa” Jacobian is such an element.) 
We show how the topology of the grid affects the complexity of the 
operators. 


1. Conservative discretization 

“Numerical methods for nonlinear conservation laws are among 
the great success stories of modern numerical analysis.” (Iserles P]) 


Conservative discretizations of partial differential equations have been 
explored for a long time. What does “conservative” mean? An early 
dehnition is due to Lax and Wendroff 0, who considered the class on 
PDEs with one spatial dimension, 

(1.1) ut +d,,{f{u)) = 0 

and called discretizations of the form 


( 1 . 2 ) 


- id: 


..., ■ ■ ■, <-j) 


At Ax 

conservative. See for an introduction to such methods. More generally, 
the formulation ( |1.1D is called conservative, and the expanded form 

Ut + f'{u)u^ 


nonconservative, with these terms carrying over to the corresponding dis¬ 
crete forms. The PDE ( p..l|) reflects, amongst other things, conservation 
of the integral of u (e.g. total mass, momentum, etc.), and the conser¬ 
vative discretization ( 0 ) preserves a discrete analog: 
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The full consequences for the discrete scheme of the form (|1.2|) remain 
unclear. 

More recently the term has been applied to PDEs that can be written 
purely in terms of intrinsic differential operators such as div, grad, and 
curl. A conservative spatial discretization is then one which preserves 
discrete analogues of these operators’ integral identities (e.g., Stokes’s 
theorem). In many cases these obey maximum principles and have ro¬ 
bust stability properties in difficult situations such as rough grids and 
discontinuous coefficients ||2^ . 


Schemes have also been developed for particular equations that inherit 
conserved quantities approximating those of the PDE. An early and fa¬ 
mous example is the Arakawa Jacobian |I|, a discretization of VxWy—VyWx 
which, when applied to the two-dimensional Euler fluid equations, pro¬ 
vides two conservation laws corresponding to energy and enstrophy, both 
quadratic functions. It is widely used in computational meteorology. 
There are many energy-conserving schemes for particular PDEs: Fei and 
Vasquez [Q for the sine-Gordon equation; Glassey 0 for the Zakharov 
equations; Glassey and Schaeffer [j^ for a nonlinear wave equation. The 
original presentations of all these are somewhat ad-hoc, the proof of con¬ 
servation relying on a telescoping sum. 

The Arakawa Jacobian has the extremely nice property that it can be 
applied to systems (in two space dimensions, with two variables) with any 
two integrals, not just energy and enstrophy. It was further explained and 
generalized to arbitrary grids by Salmon and Talley in [|T^ . It is this sys¬ 
tematic approach that we generalize in this paper to equations with any 
number of integrals, space dimensions, and variables. Our formulation 
includes all integral-preserving discretizations. 

Having integrals of course reduces the evolution to a smaller space, 
which, when their level sets are compact, gives the method a form or 
nonlinear stability. Often, more is true: Preservation of a discrete form of 
f udx hy the Lax-Wendroff form (|1.2| ) leads to correct shock speeds, and 
preservation of energy and enstrophy by the Arakawa Jacobian prevents 
energy cascading to small length scales 

The Euler equations, the sine-Gordon equation and so on are examples 
of Hamiltonian PDEs, which suggests that one should look for semi- or 
fully-discrete forms which preserve not only a discrete energy but also a 
discrete Hamiltonian (symplectic) structure. For systems with canonical 

O- 


But even 


Hamiltonian structure, this possibility was explored in 
before the importance of Hamiltonian PDEs was widely recognized, for 
which a watershed event was perhaps the 1983 conference [^|, it had 
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been pointed out by Morrison [^] that spatial discretizations of non- 
canonical Hamiltonian PDEs will not normally be Hamiltonian. One 
example apart, the curious ‘sine bracket’ Hamiltonian discretization of 
the Jacobian |2^, this is a difficult and essentially unsolved problem. 
Probably the right generalization of ‘Hamiltonian’ has not yet been found. 

We are thus reluctantly led to consider only energy-conserving dis¬ 
cretizations. Or perhaps we should not be reluctant: Simo et ah 


have argued and presented detailed evidence from elastodynamics that 
conserving energy leads to excellent nonlinear stability properties that 
preserving symplectic structure does not. (Essentially because symplec- 
tic schemes can only brake the fast modes, whereas energy-conserving 
schemes can also damp them.) 

In Hamiltonian systems, energy is normally seen as playing a distin¬ 
guished role. Yet there may be other conserved quantities just as impor¬ 
tant for the long-time dynamics. Some of them, the ‘Casimir’ integrals, 
can be due to the Hamiltonian structure itself. Non-Hamiltonian systems 
can also have conserved quantities. Even in the ODE example of the free 
rigid body, the relationship between schemes preserving energy and/or 
momentum and/or symplectic structure is quite complicated ra¬ 
in this paper we go some way towards uniting these different integrals 
and different points of view. Our goal is to develop a methodology for 
building spatial discretizations that preserve discrete analogues of any 
given set of integrals. It should be systematic, all-inclusive, and reproduce 
known schemes. We do this in a formulation in which the integrals appear 
explicitly; the integrals themselves can then be discretized in any way. 
The basic “hnite difference molecule” is now a completely skew-symmetric 
tensor. Symmetry plays a fundamental role, and we will see how the skew- 
symmetry of this tensor interacts with other desired symmetries of the 
scheme (e.g., translational and rotational invariance) in a nontrivial way. 


2. Hamiltonian and other PDEs with integrals 

We consider PDEs with independent spatial variables x G and de¬ 
pendent variables u{x) G M™'. We loosely call m the “number of vari¬ 
ables.” The relevant class of sufficiently smooth real-valued functionals 
of u{x) will be denoted U. A Hamiltonian PDE is specihed by a Hamil¬ 
tonian 7-f G W and a Poisson bracket { , } : U x U ^ U: 


( 2 . 3 ) 


ii = {u, H}, 
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where the Poisson bracket is bilinear, skew-symmetric 

(2.4) {T^Q} = -{Q^T}, 
and obeys the Jacobi identity 

(2.5) {T, {Q, n}} + {G, {n, + {n, g}} = o 

and the Leibniz rnle 


( 2 . 6 ) {J^,gn} = {j^,g}n + {j^,n}g 


for all JF, 7Y G U. In fact, these axioms imply the existence of a 
Hamiltonian (or ‘Poisson’) operator JL, such that the Poisson bracket 
takes the form 


(2.7) 




6u 


dx 


and 


( 2 , 8 ) 


becomes 


jn 

u = J— 
ou 


A particnlarly important example is the conservation law (O), which 
takes this form with x, u & R and 


(2.9) 


n 


F{u{x)) dx, J = dx, F' = f. 


The system ( p.3|) has fnnctional X G W as an integral if X = {X, 7Y}=0. 
Some integrals C are distingnished in that {C,F} = 0 for all X” G W; 
they are called Casimirs. The operator J = has a single Casimir, 
C = f udx, because ^{dC/du) = dxl = 0. 

The peculiarly Hamiltonian character of these PDEs is due to the Ja¬ 
cobi identity satished by the Poisson bracket. If, as argued previously, 
we discard this identity, what class of systems result? Energy is still con¬ 
served, because of the skew-symmetry of the bracket. Systems ( p.3|) may 
still have integrals and the operator J' may still have Casimirs. Do such 
systems retain any other special properties? The answer is no: all systems 
with an integral can be written in the form ( |2.3| ) for some choice of 
the bracket (or equivalently, for some choice of the skew-adjoint operator 
J'). So to study systems with an integral, and their discretizations, we 
may without loss of generality assume the form ( |2.3|) . 

This is most easily seen in the hnite dimensional case, as shown recently 
by Quispel and Capel [0. Take n G M"' as coordinates on phase space. 
Poisson (= non-canonical Hamiltonian) systems 


u = {u,H} = J{u)VH{u) 



SPATIAL DISCRETIZATION OF PDES WITH INTEGRALS 


5 


have integral H. But suppose an arbitrary system ii = f{u) has integral 
H. Let 2 ; = ViL and 

(2.10) Jy = 

Then J is skew-symmetric, and JVH = / as required. This J is singular 
at critical points of H, but in |]^ it is shown that if / and H are smooth, 
and the critical points of H are nondegenerate, then there is a smooth 
matrix J such that / = J'VH. 

This idea extends easily to systems with any number of integrals: The 
system of ODEs u = f{x) has integrals ... ,P if and only if there 
exists a totally skew-symmetric (p-|- l)-tensor K such that for all x where 
the vectors VT are linearly independent, 

dl^ dP 

( 2 . 11 ) U = K,^k...T^T^.... 

OXj uXk 

For, suppose K exists. Then P = f ■ VP = 0, so each P is an integral. 
Conversely, suppose / has integrals P. Then (using exterior algebra, see 


K = 


f AVP---AVP 
det{VP ■ VP) 


satishes (|2.11|) . K is determined uniquely only in the case n = p + 1] see 
131 for further details. We write the inner product ( p.ll|) as 

f = K{Vl\VP,...). 

What about Casimirs? Suppose that instead of contracting K against 
all the integrals, as in (p.ll|) , we contract against just one, say P. Then 
K = K(VP) is a skew p-tensor which has P as a Casimir, in the sense 
that K{VP) = 0 and any differential equation formed from this K (as in 
(|2.11|) ) will have P as an integral. But there are many different iC’s, all 
generating the same system ii = f, that do not have P as a Casimir. Thus 
the distinction between the Hamiltonian, other integrals, and Casimirs, 
that was present for Hamiltonian systems, is lost now. We are free to 
move between different representations of / as needed. 

(There is another importance difference. If J G satishes the 

Jacobi identity and has locally constant rank m, then J automatically 
has n — m Casimirs H- This need not be true if the Jacobi identity 
does not hold: there may not be n — m functions whose gradients span 
J’s nullspace. This is another reason for constructing J’s as above that 
automatically have the required Casimirs.) 
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One can use the tensor K to define a (p + l)-bracket, 

which is multilinear, a derivation in each argument, and completely anti¬ 
symmetric. Such brackets have been revived in modern times by Nambu 
ISl , who, amongst other things, introduced the 3-bracket on given by 
Kijk = Sijk (the alternating tensor). This gives systems of the form 

u = K{yi\ VP) = VP X VP. 

In particular, the free rigid body takes this form with P = \\u\^ being 
total angular momentum and P = \ being the kinetic energy, 

where the Ai are the body’s moments of inertia. Contracting against 'WP 
gives the standard, Lie-Poisson form of the equations, 

, , ,dP dP 

tli 


However, later studies iH 


Uk ■' ^ duk 

attempting to build a true generalization 
of Hamiltonian mechanics from such brackets, have found that not all 
constant iP’s satisfy the required “fundamental identity” (the analogue 
of the Jacobi identity); its solutions all lead to systems with n—1 integrals. 
So it is not clear that interesting dynamics as well as interesting algebraic 
structure will arise in this way. 

The situation for PDEs is formally the same. For example, for a PDF 
it = f{u) with one integral X, one can define the operator P analogously 
to ( grrUD by 

/ {f{x)z{x') — f{x')z{x)) v{x') dx' 


J{x, x')v{x') = 


f§(x')z(x')dx' 


where z(x) is any smooth function, and 

f{x) = J J{x,x')^{x')dx'. 

The extension to multiple integrals {Xfy is similar. However, questions 
of convergence of the integrals arise, and the nonuniqueness situation is 
much worse: it is not clear how to construct local operators, for example. 
However, as our goal is to construct finite-dimensional finite difference 
operators, the representation (|2.11|) is sufficient. 

Looking back at Eq. m, we see that it encompasses two impor¬ 
tant conservation laws expressed in two different ways. The fact that 
dx{Su/6u) = dxl = 0 means that the Casimir f udx is an integral, and 
the conservative scheme (|1.2|) maintains a discrete analog of this. In this 
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example, this property is relatively easy to preserve under discretization: 
a system has ^ • Ui as an integral if and only if it can be written in the 
form iii = JijFj, where J is not necessarily skew, but Jij = 0 for all j. 
Without loss of generality we can take J to be in the form (|1.2|) , with just 
two nonzero diagonals. The form of the Fj chosen in (|1.2|) is necessary 


for translation invariance. 

Secondly, if / is a variational derivative, / = 5H/5u say, then the 
skew-adjointness oi J = dx means that Ti is an integral. To preserve this 
property under discretization means taking (again, without loss of gener¬ 
ality) iii = Jij^i where J is skew symmetric. Note that such a J need 
not a priori have ^ Uj as a Casimir, and, similarly, the nonsymmetric J 
used in ( [I.2|) does not preserve any discrete H. Thus, the two expressions 
of conservation laws are in fact independent. 

In this paper we generalize the second form. The first form is decep¬ 
tively simple in this example, because the Casimir is so simple. It is not 
clear how to modify ( [I.2| ) to incorporate different Casimirs. In the second 
form the integral appears explicitly and, once J is found, any quantity 
can be conserved. 

Before continuing, we mention one trivial but complete solution to the 
whole problem. Why not contract K against all the integrals and have 
simply iii = fi, with all F being integrals of /? That is, / ■ VP = 0 
for all j. This case is already included in the above formulation, with / 
regarded as a “skew 1-tensor.” The equations / ■ VP = 0 are linear and 
can be solved in many ways, for example, by starting with an arbitrary / 
and projecting to the subspace {/ ■ VP = 0 j = 1,... ,p}. One objection 
is that this solution is so general that it is not clear how to proceed in 
any particular case. For example, to modify / as little as possible one 
might choose orthogonal projection, but this will couples all of the /*. By 
incorporating more of the known structure of the problem we can work 
more systematically. 


3. Method of discretization 


We wish to construct discretizations of the form ( p.ll| ). There are two 
ways to proceed. One could take a particular PDF, write it in the form 
ii = IC{6F ,...), and discretize this skew-adjoint operator 1C, preserving 
skew symmetry. This is difficult, if only because such formulations of 
PDFs are new and have not been widely developed yet. Instead, we study 
systems of the form ( TO in their own right, constructing elementary 
tensors Ppfc..., and seeing what PDFs they can be used to approximate. 
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That is, we establish (in a sense defined below) 

K{v ^,..., yP)... = lC{v\ ...,vP) + C>(h"). 

We call K a finite difference tensor and K{v^,, v^) a finite difference. 
Then, the integrals X* can be discretized in any way, say by 

r{u) =lfiu) + 0{h'^) 

(this amounts to a numerical quadrature) giving the conservative system 
of ODEs 


(3.12) 




Since this form includes all systems with integrals /*, we can be confident 
of not missing any in our construction. 

We start with an elementary example illustrating how easy it is to 
break skew symmetry. With one integral in one space dimension, we are 
seeking an antisymmetric matrix K. On a constant-spaced grid, central 
differences have 


(3.13) 


K = ^ 

2h 


1 

0 


which is antisymmetric. On a non-constant-spaced grid, if we let p{x) 
be the quadratic interpolating {xo,vo), (xi,ni), and {x 2 ,V 2 ) and use the 
estimate v[ = fi{xi), the associated matrix K is not antisymmetric—it 
even has a nonzero diagonal. Taking its antisymmetric part is not a good 
idea, as we have no idea what operator approximates. Indeed, it is 
not immediately clear what bandwidth is required to achieve order 2, say, 
with an antisymmetric matrix. The element which is relevant to 

■Uj, must also contribute to iii-i. 

Nonconstant operators also pose problems. Let /C = udx + dxU. Central 
differences are not skew, but it is not obvious that the skew matrix 


(3.14) 


/ 


K = 


\ 


'^i—l 


0 Ui 

-Ui 0 Ui+i 


•• / 


is a discretization of /C, or how to increase its order from 1. 
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Below we develop some requirements on the tensors K, and construct 
all the elementary ones, for various numbers of integrals, space dimen¬ 
sions, and grids. 


4. Definitions & theory 

The fundamental objects are the grid L, the index set M, the symmetry 
group G, and the skew tensors K G which we now dehne. 

Let L be a set of indices of grid points. To each index i G L there 
is a physical point Xi G Let {1,... ,m} be the set of indices of the 
dependent variables, so that the full, discrete state space is indexed by 
the index set 

M = L X {1,..., m}. 

A grid function is a real function on M; for example, the system state is 
given by the grid function u : M —> M. Its value at point z = (/, a) G M 
is written Ui = ug^a)- (That is, we are assembling all the unknowns into 
a big “column vector.”) When m = 1, we drop the second subscript 
entirely. 

For simplicity, we only consider the interpretation of this function in 
which ~ Ua{xi). (Staggered grids and Ui representing other func¬ 
tionals of u{x) do not affect our main line of argument.) 

The p functions, and their corresponding grid functions, which are to 
be inserted in ( p.l2| ) are denoted ... ,v^. However, when p = 1 we 
denote it u, and when p = 2 we denote them v and w, to reduce the 
number of indices. 

A discretization of a PDF is thus a vector held on M^, and a discretiza¬ 
tion of an operator /C is a skew (p -|- l)-tensor KiQ...ip = Ki = ^ 

ij G M] i.e., K G A^''‘^(M^). Since we are trying to construct such 
tensors, intermediate steps will also involve nonskew tensors, i.e., real 
functions on 

K approximates JC to order r if 

.. .u^(xip) = lC(ui,... ,uP)(xiJ -F 0{W') 

for all smooth We sometimes drastically abbreviate this to 

Kv = )Cv + 0{lf). We also abbreviate vl_^ • • • pfj, to P-i- 

Let the grid L have a nonnegative distance function \j — k\. This 
extends to M by |(j, a) — (/^,/d)| := |j — k\. The bandwidth of K is the 
smallest c such that = 0 for all i such that \ij — ik\ > c. 

Two examples are the Euclidean distance \xi — Xj|, giving the “Eu¬ 
clidean bandwidth,” and the minimum number of edges traversed going 
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from i and j, where the grid points have been connected to form a graph, 
giving the “graph bandwidth.” For example, K in Eq. ( p.l4[) has graph 
bandwidth 1 . iF is local if it has a hnite bandwidth even on inhnite grids. 

The grid L may be structured, like a square or triangular lattice, or 
unstructured. Let G be a symmetry group acting onM = Lx{l,..., m}. 
We usually consider only spatial symmetries, those which are the identity 
in their second slot, i.e. 7^29(1 ,«) = These merely rearrange grid 
functions on L. 

Furthermore, G is usually a subgroup of the symmetry group of the 
continuous physical space. For example, suppose this space is the plane. 
Many PDEs of physical interest are invariant under the group E{2) of 
Euclidean motions of the plane (the semi-direct product of rotations, re¬ 
flections, and translations, 0(2) (g)M^). A discrete version of such a PDE 
can inherit some of this invariance if L has a subgroup of E(2) as a sym¬ 
metry group. Examples are the square lattice, which has O 4 (s) as a 
symmetry group (8 rotations and reflections, plus discrete translations) 
and the equilateral triangular lattice, which has O 3 ©Z^. For the Eu¬ 
clidean group of the sphere there are no such natural lattices, and the 
dislocations that occur, e.g., when triangulating an icosahedral grid, are 
known to destabilize numerical methods and create artifacts in the solu¬ 
tions. 

Some nonspatial symmetries can also be included. For example, for 
a PDE involving f{v 2 x — viy) we might include the map {x,y,vi,V 2 ) ^ 
{y,x,V 2 ,vi) in G. 

Order of accuracy can also sometimes be expressed as a symmetry. 
One way to ensure second-order accuracy is for the expansion of the 
discretization error in powers of the spatial grid size h to have only odd 
or only even terms present. This is equivalent to being invariant under 
the operation h 1 —*• —h, or x 1 —*• —x. This can only apply if x 1 —*• —x 
is a symmetry of the lattice itself, which it is for square and triangular 
lattices. 

To include this possibility we equip each element of G with a sign, 
sgn( 5 f) = ± 1 , such that G is homomorphic to Z 2 . The map corresponding 
to h I— >■ —h would then have sign 1 (—1) when the operator has an even 
(odd) number of derivatives. 

The action of G extends to an action on AP’''^(M^) by 


which we write as 


{gK)i = g{Ki) = 
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A tensor K is G-invariant if gK = {sgng)K for all g E G. 

Thus we have the following requirements on the hnite difference oper¬ 
ator K: 

• K should be completely skew-symmetric; 

• K should be G-invariant; 

• K should be as simple as possible; 

• K should approximate the desired continuous operator to the desired 
order; 

• K should be local. 

However, these requirements conflict with each other. 

One way to construct operators such as K with the required symmetry 
properties is to sum over the symmetry group. Given any tensor K, 

(4.15) ^ sgn(a) sgn(^) JC(g(i)) 

is completely skew-symmetric and G-invariant. This suggests two ways 
to construct symmetric K's\ 

• Start with a K which approximates the desired continuous operator, 
and symmetrize it; 

• Start with a very simple K, such as a basis element for the space of 
[p + l)-tensors, symmetrize it, and see what continuous operator it 
approximates. 

A major drawback of the hrst strategy is that we have no control over 
what the symmetrized K approximates. 

The second strategy builds a “library” of all such difference operators, 
from which linear combinations can be taken as desired. However, the 
form (|4.15 ) is not convenient for writing down these operators in the 
usual way, which requires the coefficients of each v]. appearing in the 
resulting grid function at a particular point That is, we want to know 
for a particular We derive hnite differences in this form in 
three stages: hrstly, for m = 1 variable; secondly, for m > 1 variables 
with no unknowns at the same grid point coupled; thirdly, the general 
case, m > 1 variables with arbitrary coupling. 

Case 1. m = 1 variable. When m = 1 we drop the second component 
of the index i = [I, a). We can take io = 0. Fix a multi-index i G 
where io = 0 and start with the elementary tensor dehned hy Ki = 1 , 
Kj = 0 for all j ^ i. We assume that the indices in i are distinct, for 
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otherwise skew-symmetrizing K would lead to the zero tensor. Skew- 
symmetrizing K gives a tensor of bandwidth maxj^k \ij — ik\- It is 

Hi= ^ sgn(p) sgVL{h)Kh{p{^i)) 

pGS'p-i-i ,/iGG 

so the vector held at the point 0 is where Iq = 0. Since K is a 

discrete delta function, there is only one nonzero term in this sum, i.e., 

(4.16) ^Hivi = ^sgn(p) sgn{h)vi, 

I p,h^l 

where the sum is taken over all p, h, and I such that p{h{l)) = i and 
lo = 0. Therefore, I = p~^{h~^{i)). Let g = h~^, a = p~^ so that Eq. 
(|4.16|) becomes 

U17) E sgn(a) sgn(^)u<,(g(i)) 

o-,g 

where the sum is over all a and g such that a{g{i))o = 0. 

Since no two indices in i are the same, for each such g^ let a be such 
that a{g{i))o = 0; then the remaining a’s that satisfy this equation he 
in S'p, the permutations of the last p indices. The sum over Sp can be 
evaluated to give a determinant, giving the vector held at point 0, 

(4.18) F{i) := sgn(a) sgn(^) det I4(g(i)) 
where 

(4.19) Gi = G G : 3cr G S'p+i such that (T{g{i))Q = 0}, 
and the p x p matrix V has (j, k) entry 

(4.20) {Vi)jk = vl, l<J,k< p. 

(The hrst subscript of /q = 0 and does not appear in the matrix.) 

We introduce a graphical notation for formulas such as Eq. (|4.18|) . 
An arrow connecting gridpoints ii,... ,ip will indicate a term det Vi. For 
p > 1, the sign factors can be incorporated by applying a permutation 
of sign sgn(cr) sgn(p) to the ij (for p = 2 and p = 3 we merely change 
the direction of the arrow if the required sign is —1, equivalent to writing 
the columns of V in reverse order) or by choosing a in ([4.191) so that 
sgn(cr) sgn(p) = 1. The reader is encouraged to refer immediately to Fig. 
([^(a)) and its associated hnite difference Eq. o to see how easy this 
is. 
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Thus, constructing skew finite differences amounts to choosing an ini¬ 
tial arrow, finding the group Gi, and finding the image of the initial arrow 
under G*. 


Case 2. m > 1 variables, distinct points coupled. Let i = {I, ex) is 

the chosen basis element, where Iq = 0. “Distinct points coupled” means 
that Ij 7 ^ Ik for all j, k. 

The single element i will contribute to the vector field at (0, aj) for 
all j] therefore, we will construct, not a single basis vector field, but the 
family of vector fields spanned by K^i a) for all cx. That is, we allow 
coupling of all components right from the start. We introduce the rank 
p-|-l, dimension m tensor G and start with Ki ^ = Th, Kj at = 0 

for j 7 ^ i. 

Passing from Eq. ( [4.171 ) to Eq. (|4.18|) only required that the Ij be dis¬ 
tinct; therefore the symmetrized vector field at the point 0 (corresponding 
to Eq. (|4.18|) in the single variable case) is 

d 

(4.21) ^ sgn(a) sgn{g)Ta det V;(g(i)) 


geGi 

ao,...,ap 


du, 




We want to find the ao-component of this vector field. To do this we 
relabel the dummy indices cx by applying in the second slot only 

to get the vector field at point 0 in component ao, 


(4.22) F{i) : = E Sgn(cr) sgn(5()Tg-i(^-i(„)) det V(a(g(z)),«) 

ai,...,ap 


Notice that in ([4.22|) , each determinant involves the same components 
of the u*. Also, if is a spatial symmetry, then 5 f“^((T”^(Q:)) = cr“^(Q:). 

The diagram notation extends easily to (|4.22|) . To the arrow a{g{l)) 
we attach the label g~^{a~^{(x)) indicating the T-tensor attached to that 
determinant. 


Case 3. m > 1 variables, arbitrary coupling. Equality of some of 
the Ij affects the sum over permutations in Eq. (|4.17|) . Let j = h{p{l)) 
where jo = 0- Let n{j) be the number of O’s in j. Then the subgroup of 
Sj+i leaving cr(j)o = 0 is not Sp as it was before. It is convenient to have 
a sum of determinants of all p x p matrices, so we express this subgroup as 
the product of Sp and the flips (0/c), k = 1,... ,n{j) — 1. This subgroups 
of Sp+i contains all permutations of the first n elements; summing over 
these merely skew-symmetrizes T. We could have imposed this in the 
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first place for simplicity. Since there will usually be Qj (a translation, 
say) with gj{lj) = 0, this is true for any set of equal elements of 1. 

Let Til = {a : <7(0 = 0 be the symmetry group of Z. To sum up, we 

(4.23) take T to skew-symmetric under Ti. 


With this assumption, each flip (0/c) gives an equal contribution and we 
can evaluate the sum over permutations to give the symmetrized vector 
held at the point 0 

(4.24) F{i) : = E sgn(a) sgn{g)n{g{l))Tg-^^-^^)) det V(^(g(j)),«) 

g&Gi 

ai,...,ap 


In the diagrams, to the arrow a{g{l)) we attach the weight n{g{l)). 

To summarize, the hnal hnite difference evaluated at the point (0, cio)) 
is given by Equations ( |4.24|) ( |4.19| ), and ( |4.20|) . Eq. (|4.24|) specializes to 
Eq. ([4.22|) when all elements of Z are distinct, and specializes further to 
(|4.18|) when m = 1. In practice, from a diagram one writes down the 
hnite diherence directly from its diagram. 

We develop these diagrams and study the resulting diherences for dif¬ 
ferent numbers of integrals and dimensions of phase space, and diherent 
symmetry groups G. 


5. Examples 

Case 1. p = 1 integral, cZ = 1 space dimension. By scaling it is 
sufficient to consider i = (0,1). The bandwidth is 1. The group Gi 
has the single element z i—^ i — 1. The permutation which brings 0 to 
the front is (0, —1) i—*• (—1,0), with sign —1. Thus we get the standard 
central diherence 

F(0,1) =vi -n_i. 

On a grid with constant spacing h, 

vi — V-i = 2hvx + 0{h^) 


All other examples are related to this one: 

1. by Richardson extrapolation, 

8F(0,1) - F(0,2) = 12hvx + 0{h^). 

2. eliminating the leading order term(s) gives hnite diherences approx¬ 
imating higher-order diherential operators: 

F(0, 2) - 2F(0,1) = 2hhxx. + 0{h^). 
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3. taking a linear combination of these basis elements gives stencils 
that approximate other hrst-order differential operators. This is 
equivalent to multiplying by a symmetric tensor To get a smooth 
continuous limit we can take, e.g., Sij = q{xi,Xj,Ui,Uj), where q is 
symmetric in its hrst and second pairs of arguments. If K is the 
tensor corresponding to the vector held F{0, 1), 

(5.25) SijKijVj = 2h{sdx + dxs)v + 0{h^), 


where s := q{x,x,u,u) (no sum on i). 

We consider this last example in more detail, since it gives the class of 
all skew tridiagonal (i.e., bandwidth 1) hnite differences. 

Firstly, suppose we have the non-constant operator J' = ud^ + d^u. 
Eq. (|5.25|) discretizes this if s{x,u) = u. The only constraint is the the 
tensor must be symmetric, to maintain the overall skew-symmetry of 
the hnite diherence. For example, we can take Sij = {ut + Uj)/2 (only 
elements with \i—j\ = 1 are actually used). This gives the hnite diherence 
tensor 


(5.26) 

K 


( 


\ 


1 

2 


-Mo - Ui 


0 

-Ml - M2 


Ml M2 

0 U2 + U^ 

-M2 - Ms 0 Ms M4 


V 


/ 


showing how the skew-symmetry is maintained. Q 

Secondly, suppose we wish to diherence on an irregular grid, where the 
data are known at the points c{xi) = c{ih )—we know v{c{xi)). Then we 
want an approximation of v' = Vc- From the chain rule, this is equal to 
Vx/cx- We cannot get this by applying (|5.25|) to m, since the terms in v 
only cannot cancel. 

We apply ( b.25|) to a function w{x,v{c{x)). This gives the equation 
{sdx + dxs)w{x,v{c{x))) = {V^dxV^)w{x,v{c{x)) = 

c[x) 

with solution 


s = 


2 c' 2 ’ 


w = cv. 


^Interestingly, the choice Sij = y^UiUj actually gives a Poisson K. This is because 
it is the image of the standard central difference under a change of variables which 
sends dx into udx + dxU. 
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For the discretization, any Wj and any symmetric Sij with this con- 
tinnous limit can be taken. If c' is known analytically, we can use the 
midpoints of the intervals to get 

1 -1 


^ 0,1 — 


2 (c'J 


2 ’ 


50,-l — 


2 (c'_J 


2 ’ 


Vn = 


c'lVi 
J 2 


C iV-i 


0{h^ 


4h i c'l c'_ 

2 2 

—that is, a second order, 2-point, anti-symmetric discretization of the 
derivative on smooth grids. With constant spacing, c{x) = x and it 
reduces to the standard central difference. 

If d is not known, we can approximate it symmetrically with 

-2 


foi — 


C-i Cq 

h 


{Vi + Vi+i){Ci+i - Ci) + {Vi + - Ci_i) 


which has bandwidth 2. It is not consistent on rough grids, however; the 
error is 0{hdd"') = 0{h~^) near a discontinuity of d. 

With m > 1 components and bandwidth 1, the two possibilities are 
/ = (0,0) and I = (0,1). In the hrst case, the tensor T must be skew 
symmetric in its only two slots, but there are no group symmetries. In 
the second case, Gi has two elements, the identity and left translation. 
Under left translation, (0,1) i—>■ (—1,0); applying a = (01) of sign —1 
maps (—1, 0) I—>■ (0, —1). In the last step we apply the a~^ to the indices 
of T. Combining both possibilities gives the hnite difference 

-T^n_i + Jvo + Tvi = {J + T- T^)v + h{T + + 0{h^), 

where J = —J"^. Further imposing the symmetry i —^ —f, of sign 
— 1, is equivalent to taking T = T"*"; then the hnite difference is sec¬ 
ond order. As above, nonconstant operators are approximated by taking 
Ti = T{xi,xi+i,ui,ui-^i), symmetric in its hrst and second pairs of argu¬ 
ments, and Ji = J{xi-i,xi,xi-^i,ui-i,ui,ui-^.i), symmetric under (13) and 
(46). 


Case 2. p = 2 integrals, d = 1 space dimension. With m = 1 
component, the simplest tensor has base index i = (0,1,2). This will 
lead to bandwidth 2. We start with the arrow 1 —> 2 (see Figure p. 
Translating left by 1 and rotating indices right (an even permutation) (i.e.. 
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-2 


-2 


-1 

- • 


bca 

abc 


(a) 


cab 

cab 

(b) 


abc 

bca 


Figure 1. Case 2, Two integrals, one space dimension, 
(a) One variable; (b) m variables. See Eq. (p.27|) for the 
finite difference interpretation of (a). 


bca 


G 


abc(2) 

(a) 


abc(2) 


o 


bca 


(b) 


Figure 2. Case 2, Two integrals, one space dimension, m 
variables, I = ( 0 , 0 , 1 ). 


performing ( 0 , 1 , 2 ) i—^ (— 1 , 0 , 1 ) i—*• ( 0 , 1 ,—!)) gives the arrow 1 —*• — 1 . 
Repeating gives the arrow —2 —*• —1. The resulting diagram is already 
symmetric under h —> —h, so we do not need to add this operation. The 
diagram in Fig. |^(a) corresponds to the finite difference 

(5.27) 

F( 0 , 1 , 2 ) = det Fi ,2 + det Vi,-i + det R- 2,-1 


= {viW2 - V2W1) + {viW-i - V-iWi) + {V-2W-1 - V-1W-2) 
Expanding in Taylor series, this is 

h^{3{v'w" — w'v") + 2{vw"' — wv"')) + 0{h^) 


With m > 1 component there is essentially one finite difference each 
with bandwidth 0, 1, and 2. We write = Tabc- 

With I = (0,0,0) (bandwidth 0), Eq. ( 4.’23|) says we must have Tabc 
completely skew-symmetric. One might not call this a “difference,” since 
it only acts on vq and tco- 

With I = (0, 0,1) (bandwidth 1), Eq. ( |4.23| ) says we have Tabc = —Tbac- 
The weight of / is n{l) = 2, since it has 2 zeros. Gi has two elements, 
the identity and a left translation. Applying the left translation followed 
by a shift-right permutation a (of sign 1 ), ( 0 , 0 , 1 ) 1 —>■ (— 1 ,— 1 , 0 ) 1 —>• 
(0, — 1 ,— 1 ), giving the arrow —1 —>• —1 with label cr“^(a 6 c) = bca. 
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(0, —1, —1) has one zero, so the weight of this arrow is 1. Together we 
get two arrows, with diagram Fig. §(a) and hnite difference 

= {Tabc + Tbca + Tcab)VbWc + h{Tabc “ Tbca)Vhw'^ + h{Tcab “ TbcaW^Wc + C>(/l^). 

(The initial, skew-symmetric term conld be removed by a term T(0, 0, 0).) 

The resulting tensor is not invariant under g : i —i. Its image under 
g is shown in Fig. |^(b). These two diagrams can be added or subtracted 
to get a tensor that is invariant with sign 1 or — 1 , as desired. 

With I = (0,1,2) (bandwidth 2), Eq. ([4.23|) says that T is arbitrary. 

The group Gi has three elements: the identity, and a shift left by 1 or 2. 
Apply the two translations gives the diagram Fig. |l](b). However, unlike 
this case with m = 1 , this is not invariant under i i—^ —i, i.e., it does 
not give a second-order hnite difference. Applying this symmetry gives 
the second row of labels in Fig. |5(b). (For example, under (0,1,2) i—^ 

( 0 , — 1 , — 2 ) the arrow 1 —> 2{abc) maps to the arrow —1 —*• —2{abc) with 
sign —1, or —2 —> —l(afec) with sign 1 . Adding these makes Tabc = Tbca, 
i.e., we can take T to be symmetric under even permutations. 

Case 3. p = 1 integral, d = 2 space dimensions. With p free indices 
in K we can only couple unknowns which span a p-dimensional subspace 
of M'^. This is equivalent to the case d = p. For example, on a square 
grid in F(( 0 , 0 ), ( 0 , 1 )) = hvy -\- 0{h^). 

Thus, to get fundamentally new hnite diherence tensors, we need p > d. 

Case 4. p = 2 integrals, d = 2 space dimensions. Consider m = 

1 and a square grid. The simplest index set i we can take is i = 

((0, 0), (1, 0), (0,1)), as shown in Fig. H(a). Unfortunately, this has lattice 
bandwidth 2 and Euclidean bandwidth a/2, an unavoidable property of 
the lattice. 

Applying the two translations gives the diagram ^(b): the simplest 
translation-invariant skew tensor. It gives h{vxWy — VyW^) + 0{h‘^). That 
is, it is an “Arakawa”-type Jacobian. 

With even operators on this grid, rehections (necessary for second or¬ 
der accuracy) and rotations coincide, which reduces the complexity of 
the generated hnite diherence. Applying them gives the diagram ^(c), a 
second-order Jacobian. Finally, applying the rotations by 7 r /2 gives the 
diagram ^(d), a Jacobian with the full symmetry group (s) T?. As 
can be shown by expanding the entire hnite diherence. Fig. ^(d) is the 
Arakawa Jacobian (hrst derived in [Q.) 
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Figure 3. Case 4, Step-by-step construction of the 
Arakawa Jacobian. 


(0,1) 

• • 

• • 

• - 

• 

bca 

• - • 

\ (1,0) 

/ \ 

/ 

\ 

ca^ \^abc 

• • • 

• • • 

• • 

• 

• • • 

• • 

•-• 

\ 

•- 

/ 

• 

abc'\ /cab 

•-• 

(a) 

(b) 

(c) 


bca 

(d) 

Figure 4. 

Case 4, The Arakawa Jacobian 

on a 

triangular grid. 


We could have stopped at Fig. ^(c); its anisotropy may be irrelevant for 
some problems, and its complexity is half that of ^(d )— 12 terms instead 
of 24. 

Consider the same problem on a regular triangular grid. Now i = 
((0,0), (1,0), (0,1)) (Fig. §(a)) will give a graph bandwidth of 1, not 2. 
Applying the two translations gives Fig. ^(b), and reflections Fig. ^(c), 
which already has the full symmetry of the grid. Thus Arakawa-type 
Jacobians are naturally suited to triangular grids. (Notice that Figs. 
§(c) and 11 (c) are essentially the same.) 

There are two points to learn from this: 
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( 0 , 1 ) 



( 1 , 0 ) 


Figure 5. Case 5, 3 integrals in 2 dimensions 


1. With p integrals, grids with p + 1 mntnal nearest neighbours around 
a cell will give tensors of bandwidth 1. This is only possible in 
dimension d > p. 

2. On some grids, the (optional) spatial symmetries coincide with some 
of the (required) skew symmetries and/or reflection symmetry (needed 
for second order accuracy). 

Cases 5 and 6 illustrate these points. 

(An m-variable analogue of the Arakawa Jacobian is shown in Fig. 
^(c). It approximates a complicated second order operator, but if Tabc is 
symmetric under even permutations, it is Sy/Sh'^TabcJ^{vb, Wc) + 0{h^).) 

Cases 5 & 6: p = 3 integrals in 2 & 3 dimensions. The above ob¬ 
servation suggests that in two dimensions, the square grid, with 4 vertices 
around each cell, is better suited to the case of three rather than two in¬ 
tegrals. With i = ((0, 0), (1, 0), (1,1), (0,1)), applying the 3 translations 
only gives a tensor which is ZJ 4 -symmetric (Fig. |^). It equals 

4:h\v^J{v^, v^) + v^J{v\ v^) + v^J{v\ v^)) + 0{h^) 

where J is the Jacobian. Taking = 1, for example, recovers the 
Arakawa Jacobian Fig. 0(d), and shows that the Arakawa Jacobian also 
has the (Casimir) integral 

It also suggests that in three dimensions with three integrals, a face- 
centered-cubic grid (the red points in a red-black colouring of a cubic 
grid) is suitable. Each vertex is surrounded by 8 tetrahedra. Tak¬ 
ing I = ((0, 0, 0), (0,1,1), (1, 0,1), (1,1, 0)) (i.e, coupling the unknowns 
around one of the tetrahedra) leads to a fully symmetric discretization 
of the three-dimensional Jacobian /dxj). Using a cubic grid with 

I = ((0, 0, 0), (0, 0,1), (0,1, 0), (1, 0, 0)) leads to a 3D Jacobian with twice 
the complexity. 
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6. Discussion 


We have presented a systematic method for discretizing PDEs with a 
known list of integrals. Since all vector helds /j with integrals ... ,P 
can be written in the form (|2.11|) , the vector fields F{i) span all integral¬ 
preserving discretizations. The required symmetry properties of K make 
the the finite differences unavoidably complicated, but sometimes the (op¬ 
tional) spatial symmetries G coincide with the (compulsory) skew sym¬ 
metries Spj^i, reducing the overall complexity of the finite difference. 

We close with some comments on future directions. 

1. We have not yet mentioned time integration. It may not as crucial 
to preserve integrals in time as in space; this is not usually done 
with the Arakawa Jacobian, for example. If it is important, we note 
that linear integrals are preserved by an consistent linearly covariant 


method (such as the Euler method used in ( |1.2| )); quadratic integrals 
are preserved by some Runge-Kutta methods such as the midpoint 
rule; and any number of arbitrary integrals can be preserved by a 
discrete-time analogue of p.l2|) [p!3| . With one integral, a simple 
method is based on splitting K |T^ . 

2. To get simpler finite differences, some of the spatial symmetries can 
be broken, as for example in the half- and quarter-size Arakawa Ja- 
cobians in Fig. ^(b,c). How important is this in practice? These 
Jacobians are still fully translation invariant. Breaking this symme¬ 
try gives even simpler tensors K, in which, e.g., different differences 
are applied to red and to black points. Is this useful? 

3. Such broken symmetries may be partially repaired “on the fly” dur¬ 
ing the time integration. At the nth time step we use the finite 
difference tensor sgn{gn)gK, with gn ranging over the symmetries. 
This decreases the symmetry errors by one power of the time step 
P, which, with At = (Ax)”, may be plenty. Most drastically, K 
could be only first-order accurate, improving to second through the 
time integration. This would require a careful stability analysis. 

4. Although our discretizations are not Hamiltonian (unless /C is con¬ 
stant), they can be volume-preserving. The system (p.ll|) is vol¬ 
ume preserving for all P if '^^dKij^ /dui = 0 for all j, k,.... 
In simple cases this is simple to arrange; K in Eq. ( 3.14 ) is not 
volume-preserving, but K in Eq. (|5.26D is. Incorporating volume- 
preservation in general is more difficult; see for a discussion. 

5. We have deliberately avoid mentioning boundaries and the precise 
degree of smoothness required of the arguments that make /C skew. 
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These are studied in [^]. If the PDE develops shocks, a careful 
weighting of the F{i) will be required to capture them well, the 
analogue of the many methods for choosing H in (p..2|) [^. The 
present work applies to inhnite or (trivially) to periodic domains. 
With hnite domains, one can start with Ki at an interior point 
i, and extend it to the boundary by skew-symmetry, giving hnite 
difference tensors satisfying certain “natural” boundary conditions. 
Also on hnite domains, there is the possibility of using global (e.g. 
spectral) methods. Of course, these are in the span of our basis, but 
that is not the best way to view them. Preserving integrals with 


global methods is studied in |14 


6 . We have concentrated on constructing skew-symmetric tensors ap¬ 
proximating skew operators. Exactly the same technique can be 
used to construct symmetric tensors. We replace the canonical sign 
function on Sp+i by any sign function a that makes S'p+i = Z 2 . If 
a{i) = 1 for all i, for example, the resulting K is completely sym¬ 
metric, and when contracted against any p — 1 of the P, has real 
eigenvalues. If negative dehnite, the P decrease in time. What is 
the relationship with the support operator method ||20||? 
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Appendix 


Consider the operator K, = d^- We were puzzled by the following: the 
tensor K in Eq. ( b.l3| ), 

/ •• \ 


2h 


0 1 

-1 0 1 


V / 

preserves not just the integral H it operates on, but also C = 
because C* is a Casiniir of K. But this two-integral discretization does 
not arise from any of the rank 3 skew-tensors we derived in Section 5, 
Case 2—Eq. ( 5.27 ) in particular. Contracting with the required integral 
C gives a discretization of d^xxj not of d^- The same is true for any other 
basis element. 

To force C to appear explicitly in the discretization, we first find a 
skew differential operator J'{u,v,w) such that J'{u,v,SC) = IC{u,v) for 
all Ti. If we restrict to a finite domain D so that 1 dx is finite, a natural 
solution is 

J (u, V, w) = uVx / w dx — uWx / V dx + u / vw^ dx. 

Jd Jd Jd 

This is a non-local differential operator, which is the resolution of the 
paradox. It only reduces to a local operator when tc = 1. Discretizing its 
derivatives by central derivatives, and integrals J wdx hy gives a 

non-local skew 3-tensor J such that J(VC) = K. In Section 5 we only 
looked at local tensors. 

It seems unlikely that the telescoping sum which makes this example 
work will work for nonlinear Casimirs. On the other hand, it is quite 
hard to destroy linear ones. Therefore we suggest the following strategy: 
temporarily disregard any known linear integrals (mass, momentum etc.). 
Construct a skew tensor so as to preserve the desired nonlinear integrals. 
Then, check that this tensor has (some discretization of) the required 
linear integrals as Casimirs. 

The situation is analogous to preserving volume, a linear differential 
invariant discussed in Section 5, note 4. 
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